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EXECUTIVE  SUMMARY 


The  main  issue  being  addressed  is  as  follows:  Given  what  is  known  and  not  known  about 
a  rock  site,  what  can  one  say  about  the  uncertainty  in  the  ground-shock  environment  that 
can  be  expected  from  an  explosive  event?  To  be  sure,  the  less  (more)  data  we  have  on  a 
site,  the  greater  (smaller)  the  ground-shock  uncertainty.  But  does  a  quantitative,  rational 
procedure  exist  that  one  can  use  to  estimate  the  environment  uncertainty  that  corresponds 
to  the  degree  of  uncertainty  on  the  site? 

As  mentioned  in  our  proposal,  there  are  four  subsidiary  issues:  (1)  Characterization  and 
quantification  of  the  variability  in  the  mechanical  properties  of  a  rock  site  -  how  can  the 
variability  be  best  encapsulated  in  a  stochastic  model?  (2)  The  prediction  of  the  ground- 
shock  response  for  the  stochastic  prescription  -  is  the  state-of-the-art  in  stochastic  finite- 
element  technology  adequate  for  this  purpose?  (3)  Extraction  of  the  local  and  spatial 
uncertainties  in  the  predicted  ground-shock  environment  -  given  stochastic  realizations  of 
the  ground  shock  environment,  how  can  their  local  and  spatial  uncertainties  be 
encapsulated?  And,  most  importantly,  (4)  the  ultimate  goal  of  attributing  these 
uncertainties  to  the  site  uncertainties  -  how  can  a  quantitative,  site-environment 
uncertainties  transfer-function  be  established? 

Such  a  transfer  function  can  be  invaluable.  Site  databases  vary  from  application  to 
application  and  have  various  degrees  of  incompleteness;  those  commonly  encountered  in 
real-world  situations  are  usually  quite  sparse.  The  ability  to  quantify  the  environment 
uncertainty  based  on  the  uncertain  site  data  will  facilitate  field-test  planning,  structural 
vulnerability  assessment,  data  collection,  and  other  goals. 

Note  that  site  categorization  for  the  purpose  of  assessing  explosive  effects  is  much  more 
demanding  when  compared  with,  say,  seismic  profiling.  The  main  reason  is  that  the 
explosive  phenomenology  involves  high  stresses  and  high  strain  rates  and,  hence, 
mechanical  properties  of  the  rock-mass  under  these  extreme  conditions.  Yet  these 
properties  cannot  be  measured  directly  and  insitu  at  present;  they  can  only  be  inferred 
from  other  measurements.  Material  behavior  under  said  conditions  is  highly  non-linear; 
thus  the  effects  of  heterogeneities  and  spatial  variations  inherent  in  geological  materials 
become  much  more  important.  Hence,  domains  that  can  normally  be  considered  uniform 
and  homogeneous  for  seismic  purposes  do  not  behave  as  such  when  subjected  to 
explosive  effects. 

We  have  demonstrated  in  this  project:  (1)  How  heterogeneity  and  spatial  variability  of 
rock  properties  at  a  site,  as  evidenced  by  state-of-the-art  measurements  taken  at  the 
demonstration  site,  could  be  captured  as  random  fields  using  geostatistical  techniques;  (2) 
How  these  fields  of  mostly  low-stress,  in-situ  or  laboratory  properties,  which  are 
representative  of  the  kinds  of  measurements  being  taken,  could  be  used  to  infer  rock- 
mass  mechanical  properties  in  the  high-stress,  high-strain-rate  response  regimes  of 
interest;  (3)  How  the  rock-mass  property  fields  were  used  in  ground-shock  simulations  to 
generate  the  corresponding  ground-shock  fields,  employing  stochastic  finite-element 
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techniques;  and,  lastly,  (4)  How  local  variability  in  the  simulated  ground-shock 
environments  could  be  readily  elicited  from  the  stochastic  realizations. 


Local  variability  addresses  the  question:  What  uncertainty  could  one  reasonably  expect  in 
the  ground-shock  measurements  at  a  location  at  the  site  in  an  explosive  test.  The  figure 
below  compares  two  cutout  views  of  the  volume  of  interest.  They  give  the  (mean  ±  la) 
bounds  on  the  distribution  of  the  peak  pressure,  in  MPa,  where  mean  is  the  mean  value 
and  a  the  standard  deviation,  respectively,  of  the  20  realizations  obtained  by  finite- 
element  simulation  based  on  stochastic  modeling  of  the  site.  As  can  be  seen  from  the 
figure,  the  mean  and  a  values  vary  with  location.  At  a  given  location,  the  (mean  ±  la) 
bounds  measure  the  variability  that  can  be  expected  for  that  location  and,  hence,  the  local 
variability. 


(a)  mean  minus  1  sigma 


(b)  mean  plus  1  sigma 


On  the  other  hand,  we  have  yet  to  demonstrate  how  spatial  variability  exhibited  in  the 
same  simulated  environments  could  be  encapsulated.  Spatial  variability  addresses  the 
question:  How  much  variation  in  the  measurement  could  one  reasonably  expect  if  it  is 
made  in  a  slightly  different  location.  A  relation  between  spatial  site-property  variability 
on  the  one  hand  and  spatial  ground-shock  response  variability  on  the  other,  has  not  been 
established1. 


It  is  posited  that  the  very  same  geostatistical  models  used  to  capture  site-property  variability  may  be 
applied  to  summarize  environment  variability.  After  all,  realizations  of  the  environment  may  be  considered 
“response  measurements”  —  measurements  that  are  made  for  all  locations  of  the  site  as  opposed  to  finite 
locations  in  site  exploration.  Hence,  a  relation  between  the  geostatistical  models  of  the  site  and  those  of  the 
“response  measurement”,  may  be  established. 
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DEMONSTRATION  DETAILS 


The  demonstration  was  carried  out  using  the  Mitchell  Quarry  (MQ)  test  bed  that  evolved 
from  a  recent  site  exploration  project  sponsored  by  DTRA  to  assess  weapons  effects  on 
underground  facilities.  The  MO  database  is  noteworthy  not  only  for  its  diversity  and 
comprehensiveness,  which  are  unmatched,  but  also  for  the  inclusion  of  RPSS  soundings 
at  over  30  locations  of  the  site.  While  still  sparse  by  geostatistics  standards,  these 
soundings,  when  supplemented  by  other  data  in  the  database,  supported  the 
characterization  of  a  spatial  variability  model  for  the  site  -  and  constituted  a  key  step  in 
our  demonstration. 

Encapsulate  Material  Inhomogeneity  and  Spatial  Variability  Exhibited  by  the  Site 
Database 

It  is  common  knowledge  that  site  measurements  are  difficult  to  interpret  at  face  value  due 
to  inherent  material  heterogeneity,  spatial  variations  and  measurement  uncertainties.  For 
example,  the  Cp  ratio  (ratio  of  insitu  p-wavespeed  to  laboratory  p-wavespeed,  squared)  of 
samples  taken  from  the  four  boreholes  at  the  site  is  given  in  Figure  1.  Holes  C2,  C3  and 
C4,  in  particular,  are  in  close  proximity  to  one  another,  and  if  one  assumes  that  the  site 
consists  mainly  of  uniform,  horizontal  layers,  readings  at  these  boreholes  taken  at  the 
same  elevation  (i.e.,  for  the  same  geological  layer)  should  be  highly  correlated.  However, 
the  results  presented  in  Figure  2  indicate  otherwise.  In  fact,  the  correlation  among 
boreholes  C2,  C3  and  C4,  which  are  very  close  to  one  another,  is  no  better  than  the 
correlation  between  these  holes,  taken  in  turn,  and  hole  Cl,  which  is  much  farther  away 
from  the  group  (see  Figure  1). 
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Figure  1.  Cp  ratio  for  samples  in  boreholes  C1-C4,  plotted  versus  elevation  (ordinate). 
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Figure  2.  Correlation  of  Cp  ratios  taken  from  the  boreholes  Cl  through  C4;  correlation  is 
computed  for  samples  that  have  the  same  elevation. 

In  this  project,  we  found  that  the  macro-structure  of  what  appear  to  be  random,  spatial 
fluctuations  (Figures  1  and  2),  could  be  uncovered  using  geostatistical  models.  The  key 
variable  is  spatial  correlation,  i.e.,  the  correlation  of  material  properties  between  two 
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points  located  some  distance  h  apart;  note  that  these  two  points  do  not  have  to  be  at  the 
same  elevation,  as  in  Figure  2.  Spatial  correlation  is  expressed  in  the  form  of  a  variogram, 
which  is  a  plot  of  the  degree  of  variability  (lack  of  correlation)  as  a  function  of  the 
separation  distance,  usually  in  a  given,  fixed  direction  in  three-dimensional  space. 
Properties  for  points  close  together  are  expected  to  have  high  correlation  (low  variation) 
and  those  far  apart  are  expected  to  have  low  correlation  (high  variation).  The  innovation 
introduced  by  geostatistics  or  variograms  is  that  instead  of  looking  at  the  relationship 
exhibited  by  a  pair  of  points,  one  pair  at  a  time,  which  can  be  baffling  as  we  have  seen  in 
Figures  1  and  2,  a  variogram  collects  as  many  qualifying  pairs  as  exist  in  the  database, 
and  takes  the  average  of  these  pairs  to  enforce  the  relationship.  In  this  way  the  macro 
structure  of  the  underlying  dependency  is  uncovered. 

For  example,  the  variogram  for  the  same  borehole  data  shown  in  Figure  1  is  presented  in 
Figure  3.  As  is  generally  the  case,  correlation  (variation)  depends  not  only  on  the 
separation  distance,  but  also  the  orientation  of  the  connecting  line,  i.e.,  pairs  with  the 
same  separation  distance  h,  but  oriented  in  different  directions,  will  in  general  exhibit 
different  correlations.  The  results  shown  in  Figure  3  indicate  that  two  major  orientations 
distinguish  themselves:  The  horizontal  and  vertical  directions,  lending  support  to  a 
basically  horizontal  layering  site  structure.  Since  the  data  points  are  more  heavily 
populated  in  the  vertical  direction  than  in  the  horizontal  direction  (many  more  samples 
taken  along  the  bore  hole  than  there  are  boreholes),  the  vertical  component  of  the 
variogram  in  Figure  3  is  much  better  defined  than  the  horizontal  component. 
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Figure  3.  Variogram  of  Cp  ratio,  based  on  borehole  data  Cl  through  C4. 
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Two  quantities  in  a  variogram  are  of  special  interest:  The  asymptote  and  intercept.  With 
reference  to  the  vertical  component  of  the  variogram  shown  in  Figure  3,  the  asymptote  is 
reached  at  a  separation  distance  of  about  15  ft.  Therefore,  properties  for  points  that  are  15 
ft  or  more  apart  vertically  no  longer  correlate;  they  vary  independently.  For  this  reason, 
this  distance  is  called  the  effective  range.  In  this  particular  example,  the  borehole  data 
used  to  generate  the  figure  are  sampled  every  3  ft  in  elevation,  so  variogram  information 
for  separation  distances  3  ft  or  less  is  lacking.  Hence,  the  intercept  (variogram  value  for 
separation  distance  =  0)  can  only  be  estimated  by  extrapolation  of  the  curve,  as  indicated 
by  the  dash  line  in  the  figure;  it  is  approximately  0.5,  or  50%  of  the  total  variance  (the 
asymptote).  Reducing  the  sampling  distance  will  yield  more  data  to  better  support  the 
intercept,  but  if  we  take  Figure  3  at  face  value,  the  result  implies  a  rather  large  component 
of  inherent  “noise”  in  the  measurement.  For  this  reason,  the  intercept  is  also  called  the 
“nugget”  effect  -  random  impurities  in  a  mined  nugget 

For  most  practical  purposes,  the  three  parameters:  Intercept  (nugget),  asymptote  and 
effective  range,  define  the  variogram  and,  hence,  the  macro-structure  of  the  spatial 
variability  at  a  site.  It  is  helpful  to  remember  the  analogy  between  the  variogram,  for 
spatial  variability,  and  the  density  function  for  random  numbers.  While  each  random 
number  (properties  at  a  location)  appears  to  have  little  overt  relation  to  the  next 
(properties  at  a  different  location),  the  samples  (properties  at  site  locations)  are  actually 
governed  by  a  structure  or  model  called  density  function  (variogram).  While  variograms 
(density  functions)  are  usually  derived  from  data,  as  indicated  herein,  they  are  sometimes 
postulated  or  assumed  as  needed,  e.g.,  when  data  are  not  available.  The  important  point 
is,  once  postulated,  the  variogram  (density  function)  dictates  the  site  variability  (random 
number  sequence).  Or,  put  differently,  slightly  different  variograms  (density  functions) 
will  lead  to  slightly  different  site  variability  profiles  (random  sequences). 

The  definition  of  the  variogram  (macro-structure  of  site  variability)  improves  when  the 
spatial  database  grows.  In  particular,  we  found  that  the  RPSS  index  dataset,  Kp,  covering 
33  locations  (Figure  4)  gives  rise  to  a  much  better  defined  structure  than  the  four 
boreholes,  as  shown  in  Figures  5  and  Figure  6  for  the  vertical  and  horizontal  components 
of  the  variogram,  respectively.  These  empirical  results  could  be  further  consolidated  into 
a  model,  given  in  Figure  7.  This  model  then  defines  the  structure  of  spatial  variations  of 
properties  at  the  site  as  reflected  in  the  RPSS  dataset.  We  call  it  the  vrnodel,  for 
variogram  model,  to  distinguish  it  from  the  variograms  themselves  which  are  based  on 
site  data. 
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Figure  4.  The  RPSS  dataset,  with  soundings  at  33  locations,  covers  a  much  larger  domain 

than  the  four  boreholes. 
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Figure  5.  Variogram  structure  for  RPSS  index,  Xp;  vertical  component.  The  effect  of 
azimuth,  as  indicated  by  the  angles  notated  in  the  figure,  is  nominal.  Hence,  one  staicture 

for  the  vertical  component  suffices. 


Several  points  are  worth  noting  with  reference  to  the  figure.  The  vertical  and  horizontal 
components  have  very  different  effective  ranges,  viz.,  25  ft  and  250  ft,  respectively. 

These  numbers  conjure  a  site  with  horizontal  layers  that  are  fairly  uniform  (large 
horizontal  effective  ranges),  while  variations  of  properties  in  the  vertical  direction  are 
rapid  (small  vertical  effective  ranges).  A  small  nugget  value,  0.5  out  of  a  total  variance  of 
4.5,  indicates  that  the  background  noise  in  the  measurement  is  small. 

With  the  help  of  the  vmodel  in  Figure  7,  the  data  from  a  finite  number  of  measurement 
locations  such  as  the  RPSS  can  be  extended  to  the  complete  site.  This  is  done  in  two 
ways: 

•  By  kinging,  which  amounts  to  interpolating  and  extrapolating  the  measurements  to 
all  points  of  the  site  using  the  vmodel,  or 

•  By  conditional  stochastic  simulation,  which  amounts  to  generating  random 
realizations  of  the  site  according  to  the  underlying  vmodel  and  conditioned  on  the 
data  points. 
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Figure  6.  Variogram  structure  for  RPSS  index,  kp\  horizontal  component.  The  effect  of 
azimuth,  as  indicated  by  the  angles  in  the  figure,  is  also  negligible.  Hence,  one  structure 
for  the  horizontal  component  suffices.  Details  for  separation  distances  less  than  20  ft  are 
not  supported  by  data  samples  since  no  two  RPSS  sounding  locations  are  closer  than  this 

distance. 

Random  fields  generated  by  kriging  and  by  simulation  honor  the  measured  data  and  the 
vmodel  but  have  their  respective  unique  characteristics  (for  details,  see  the  many 
excellent  texts  available  including  Goovaerts,  19972 3). 

Continuing  with  the  illustration  based  on  the  RPSS  data  of  the  demonstration  site,  the 
“kriging”  results  based  on  the  data  and  the  vmodel  of  Figure  7  are  shown  in  Figure  8.  The 
kriging  estimates  show  a  high  kp  layer  around  elevation  620  ft  and  a  low  kp  layer  at 
elevation  640  ft.  If  one  equates  the  value  of  kp  with  rock  quality,  the  profile  given  in  the 
figure  can  be  interpreted  as  a  site  with  roughly  horizontal  layers,  composed  of  a  medium- 
rock  matrix  (yellow),  with  a  weak  layer  embedded  in  the  upper  portion  (cyan)  and  a 
strong  layer  embedded  in  the  lower  portion  (red)  . 


2  Goovaerts,  P.  Geostatistics  for  Natural  Resources  Evaluation,  Oxford  Uuiv.  Press,  NY,  1997. 

3  An  independent  study  by  Sandia,  conducted  and  completed  prior  to  our  effort  but  also  based  on 
geostatistics,  arrives  at  the  same  general  conclusions  about  the  site  albeit  with  a  slightly  different  variability 
structural  model. 
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Figure  7.  Variogram  models  based  on  the  RPSS  dataset,  i.e.,  the  curves  are  fitted  to 
empirical  data  shown  in  Figures  5  and  6,  respectively. 
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Figure  8.  The  technique  of  kriging  extends  A p  from  the  33  soundings  to  all  points  of  the 
site.  A  cut-away  view  of  the  site  volume  of  interest,  measuring  about  100  ft  by  300  ft  by 
80  ft  and  divided  into  3’x3’x2’  cells  (35  cells  in  x  @  3\  100  cells  in  y  @  3’,  and  40  cells 

in  z  @2’),  is  shown  here. 
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Extension  of  the  RPSS  sounding  data  to  the  complete  site  was  also  carried  out  by  means 
of  (conditional  sequential  Gaussian)  stochastic  simulation  technique,  and  the  results 
given  in  Figure  9.  The  simulation  results  are  more  spotty  (salt  and  pepper)  than  kriging 
results  as  expected  -  Kriging  is  an  interpolation/extrapolation  process  and  tends  to  smooth 
out  the  site  profile;  high  values  are  under-estimated  and  low  values  are  over-estimated,  as 
is  confirmed  by  the  respective  histograms  given  in  Figure  10.  On  the  other  hand,  because 
of  their  spottiness,  simulation  results  are  sometimes  difficult  to  interpret  although  in  the 
current  example  the  strong  (red)  and  weak  (cyan)  layers  are  clearly  visible  (Figure  10). 
But  contrast  them  with  the  more  prominent  layers  in  the  kriged  field  shown  in  Figure  9. 
Note  also  that  while  there  is  only  one  kriged  field  corresponding  to  the  vmodel  and  a 
method  of  kriging  (how  interpolation  and  extrapolation  are  to  be  done),  stochastic 
simulation  gives  rise  to  an  infinite  number  of  realizations,  all  conforming  to  the  vmodel 
and  the  conditioning  data  points. 


Figure  9.  Extending  the  RPSS  sounding  data  to  the  volume  of  interest  by  means  of 
(conditional  sequential  Gaussian)  stochastic  simulation.  The  figure  shows  one  of  many 

realizations  generated. 

It  should  be  obvious  that  the  generalization  of  measured  data  to  the  complete  field  via  the 
vmodel  (either  by  kriging  or  simulation)  depends  on  details  of  the  variogram.  In 
particular,  decreasing  the  nugget  reduces  the  spottiness  and  smoothes  out  the  field,  i.e., 
less  noise;  increasing  the  range  does  likewise.  Sensitivity  of  the  random  field  to  these  key 
parameters  of  the  variogram  has  been  investigated  but  details  will  not  be  shown  here.  As 
the  database  grows,  e.g.,  more  RPSS  soundings,  these  key  parameters  become  better 
defined  and  the  uncertainty  in  the  vmodel  and,  hence,  the  random  field,  decreases.  Hence, 
the  vmodel  captures  the  data  and  their  associated  uncertainties;  coarseness  of  data  is 
reflected  in  a  coarse  vmodel. 
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Figure  10.  Comparison  of  histograms  of  the  kriged  and  simulated  fields.  Kriged  fields  are 
smoother  (less  fluctuations;  smaller  variance) 

Generate  Fields  of  Rock-Mass  Properties 

As  is  well  known,  ground  shock  studies  require  characterizing  the  appropriate  rock-mass 
property  fields  for  the  site.  These  properties  include,  among  others,  porosity,  rock-mass 
rating,  rock-mass  bulk  modulus  and  constrained  modulus  to  the  stress  level  of  interest, 
rock-mass  unconfined  compressive  strength,  strength  as  a  function  of  pressure,  and  void 
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closure  strain.  In  addition,  properties  of  the  intact  rock  such  as  grain  density,  unconfined 
compressive  strength  and  elastic  moduli  are  needed  to  round  out  the  material 
prescription. 

We  have  demonstrated  how  random  fields  of  certain  insitu  measurements,  viz.,  the  A p 
index  from  the  RPSS  soundings,  can  be  generated  by  uncovering  the  underlying  vmodel 
and  using  the  vmodel  to  extend  the  data  to  the  full  site  via  stochastic  simulation.  Ideally, 
if  the  relevant  intact  and  rock-mass  properties  required  for  ground  shock  studies  are  also 
measured  to  the  same  extent  as  A p,  the  stochastic  realization  process  used  on  A p  can  be 
repeated  to  the  measured  rock-mass  properties  and  the  required  random  fields  generated. 
But  unfortunately,  direct  measurements  of  rock-mass  properties  are  scarce.  A  larger  body 
of  intact  properties  exists,  but  it  is  also  limited  to  what  can  be  derived  from  samples  taken 
from  the  (small  number  of)  boreholes,  and  not  from  the  (more  numerous)  RPSS 
soundings.  The  current  state-of-technology  is  that  rock-mass  properties  are  inferred, 
based  mainly  on  rock-engineering  practices:  they  are  inferred  from  the  intact  properties 
such  as  porosity  (cp)  then  degraded  based  on  estimates  of  rock  quality  such  as  rock-mass 
rating  (RMR).. 

That  being  the  case,  our  demonstration  followed  the  same  strategy  to  derive  rock-mass 
properties.  But  it  is  easy  to  see  that  other  practices  and  schemes  can  be  readily 
accommodated  by  the  method  and  that  feasibility  of  the  overall  approach  being 
demonstrated  does  not  depend  on  the  particulars  of  how  rock-mass  properties  are 
inferred.  The  important  point  to  note  is  that  as  with  any  inference  procedure,  inference 
introduces  an  associated  uncertainty  which  contributes  to  the  overall  uncertainty  of  the 
end  product,  i.e.,  the  predicted  ground  shock  environment.  This  uncertainty  must  be 
quantified  and  incorporated  into  the  end  product.  We  chose  to  side-step  this  uncertainty 
for  the  time  being,  so  that  we  may  focus  on  the  generation  of  the  random  fields  that 
enable  ground  shock  simulation. 

Accordingly,  our  procedure  consists  of  two  steps.  First,  (random)  RMR  and  cp  fields  are 
generated  from  the  laboratory  and  field  data  in  the  MQ  database,  mainly  the  Vp  and  Vs 
(compression  and  shear  wave  speed)  data.  Since  these  data  exist  in  various  degrees  of 
completeness,  and  are  largely  limited  to  samples  taken  from  the  boreholes,  their  random 
fields  are  augmented  by  leveraging  the  vmodel  of  the  more  extensive  measurements,  viz., 
the  RPSS  sounding  index  Xp  introduced  earlier.  Second,  the  desired  rock-mass  property 
fields  are  generated  by  mapping  from  RMR  and  cp  according  to  existing,  well-accepted 
relations  in  rock  engineering  (e.g.,  Hoek  et  al,  2002). 

Figure  1 1  gives  an  illustration  of  one  of  the  random  insitu-Vp  fields  obtained  in  the  study. 
Note  that  variability  exhibited  in  the  data,  as  defined  by  the  borehole  samples  (sub-figure 
b),  is  retained  in  the  generated  random  field  (sub-figure  c).  The  insitu-  Vp  field  is  mapped 
to  the  porosity  field  using  the  empirical  data  summarized  in  Figure  12,  ending  with  the 
result  shown  in  Figure  13. 
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Figure  1 1.  Simulated  insitu-Fp  field.  The  figure  shows  one  of  many  random  samples 

generated. 


Figure  12.  Mapping  between  insitu-Fp  and  porosity  for  limestones.  Whisker  data  are 
from  the  MQ  database.  Backdrop  data  are  taken  from:  Assefa,  S.,  McCann,  C.,  and  J. 
Sothcott,  Velocities  of  Compressional  Waves  and  Shear  Waves  in  Limestones, 

Geophysical  Prospecting,  2002. 


16 


UttJ  l  .  *  1  |  J-  -I  |  I  ^  .  r  I 

00  S  «  U  H  l« 


(b)  llisiouram  based  on  simulated  field 

Figure  13.  Porosity  field  corresponding  to  Figure  12. 

Similarly,  random  RMR  fields  based  on  the  square  of  the  velocity  ratio,  i.e.,  ( insitu - 
VpUnt&ci-VpY,  are  generated,  and  a  sample  field  is  presented  as  Figure  14.  Direct 
information  on  RMR  can  also  be  found  in  the  MQ  database,  and  can  be  incorporated  in 
the  simulation  process  as  conditions.  We  have  not  done  so  because  of  uncertainty  in 
interpreting  the  data  and  limited  resources  available. 


(8)  Simulated  field 


(c)  Histogram  based  on  simulated  field 


Figure  E-14.  Simulated  RMR  field  based  on  square  of  the  velocity  ratio. 


17 


In  the  second  step,  given  a  particular  rock  type,  e.g.,  limestone  (for  this  exercise,  we 
assumed  that  the  sandstone  &  shale  layers  at  MQ  were  limestone),  and  a  porosity  cp, 
generic  relations  enabled  us  to  look  up  the  corresponding  intact  properties.  Using  the 
RMR  information,  we  were  able  to  map  the  intact  properties  to  ms/tw-properties  and, 
finally,  to  rock-mass  properties  using  rock-engineering  guidelines.  Details  will  be 
skipped,  but  two  properties  are  singled  out  for  illustration  below:  The  constrained 
modulus  M,  and  the  unconfined  compressive  strength  UCS. 

The  intact  constrained  modulus  undergoes  a  two-step  reduction,  to  become  the  rock-mass 
constrained  modulus: 


RMR*M 


lab 


M 


insitu 


/3  =>  M 


rockmass 


where 


P  =  \ 


0.15  if  RMR  <  60 

- 1.125  +  0.02125* RMR  \{RMR>60 


is  a  modulus  reduction  factor  attributed  to  experimental  data  obtained  by  Stagg  & 
Zienkiewicz  (1974).  And  the  rock-mass  strength  is  related  to  the  intact  strength  by: 


r  *  UCS 

J 


lab 


Tyvnr ^  rockmass 


rs  =  Je(RMR~mi9  undisturbed 

This  relation  is  attributed  to  Hoek  &  Brown  (1980,  see  www.rocscience.com). 

=  ^(RMR-myo  disturbed 

J 


Sample  results  are  given  in  Figure  15: 
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(b)  Rockmoss  unconf,  comp,  strengtli  (MPa) 


Figure  15.  Sample  realizations  of  the  rock-mass  constrained  modulus  and  rock-mass 

unconfined  compressive  strength. 

Perform  Stochastic  Finite-Element  Simulation  of  Ground  Shock  Propagation 

Stochastic  finite-element  simulation  refers  to  the  repeated  application  of  a  deterministic 
finite-element  ground-shock  simulation  tool  to  realizations  of  the  random  field.  Each 
application  takes  one  realized  random  field,  performs  a  deterministic  ground-shock 
simulation,  and  generates  a  deterministic  ground-shock  response.  When  all  is  said  and 
done,  each  realization  of  the  random  field  will  have  a  ground-shock  field  counterpart. 
The  variation  in  the  ground-shock  response  from  realization  to  realization  defines  its 
uncertainty,  which  can  be  attributed  solely  to  the  uncertainty  in  media  properties  since  the 
energy  source  is  also  deterministic  -  it  is  the  same  for  all  simulations.  This  was  done  so 
we  could  focus  on  the  effects  of  media  properties. 

Deterministic  finite-element  simulation  of  ground-shock  propagation  is  quite 
straightforward  given  the  prevalent  use  of  finite-element  tools  in  the  community.  The 
following  may  be  note-worthy: 

Assignment  of  rock-mass  properties  to  the  elements  -  The  volume  of  interest  is 
discretized  into  three-dimensional  brick  elements,  each  of  which  is  assigned  rock-mass 
properties  that  correspond  to  a  realization  of  the  random  field  described  previously.  The 
rock-mass  properties  are  uniform  within  an  element,  but  the  rock-mass  properties  of  a 
given  element  will,  in  general,  vary  from  realization  to  realization  as  dictated  by  the 
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stochasticity  of  the  random  field.  The  element  size  is  an  important  issue,  especially  in 
relation  to  the  cell  size  used  in  generating  the  random  fields;  the  element  size  may  be 
different  from  the  cell  size.  For  simplicity,  we  chose  an  element  size  that  is  the  same  as 
the  cell  size,  viz.,  3’x3’x2\ 

Implementation  of  properties  in  a  constitutive  model  -  The  next  step  in  the 
simulation  is  to  implement  the  rock-mass  properties  of  an  element  in  a  constitutive 
model,  so  that  it  can  be  exercised  by  the  ground-shock  code.  The  choice  of  constitutive 
model  (and  its  computer  implementation)  is  usually  constrained  by  the  ground  shock 
code  used;  each  finite-element  tool  has  its  own  material  modeling  preference.  Our  use  of 
FLEX4  in  this  demonstration,  for  example,  required  that  the  rock-mass  properties  be 
implemented  in  a  SFT1  model5.  Suffice  to  say  that  the  same  set  of  rock-mass  properties 
may  be  implemented  in  different  constitutive  models  for  use  in  various  ground-shock 
codes,  and  the  various  implementations/codes  will  not  in  general  produce  the  same  result. 

For  our  purpose,  a  set  of  SFT1  parameters  was  generated  to  correspond  to  a  particular  set 
of  prescribed  rock-mass  properties.  The  SFT1  parameters  are  then  incorporated  as  one 
“material  type”  in  the  material  library,  so  that  they  can  be  exercised  in  a  simulation. 
The  determination  of  the  SFT1  parameters  from  prescribed  rock-mass  properties  is  not 
trivial;  it  is  a  source  of  uncertainty  in  its  own  right.  Furthermore,  since  rock-mass 
properties  have  continuous  valuation  within  their  respective  ranges  of  variation,  we  have, 
in  effect,  an  infinite  number  of  material  types  in  the  library.  This,  obviously,  cannot  be, 
and  we  circumvented  the  difficulty  as  described  in  the  following. 

(a)  In  order  to  keep  the  number  of  materials  reasonable,  we  picked  10  discrete  values 
of  RMR  (0.5<RMR<1.0),  and  10  discrete  values  of  porosity  (0<  cp  <0.2).  This 
resulted  in  100  unique  (RMR,cp)  pairs  and  hence  100  unique  materials.  For  each 
cell  in  each  realization,  we  simply  assigned  the  material  with  the  closest  (RMR, 
cp). 

(b)  To  construct  nonlinear  SFT1  models  from  the  (RMR,  cp)  values  we  adopted  the 
rock  engineering  correlations  recently  published  by  Chitty6.  Chitty’s  report 
became  available  during  the  course  of  this  investigation.  It  includes  many  of  the 
accepted  rock  engineering  relations  mentioned  above,  and  extends  them  into  the 
higher  pressure  regime  relevant  for  shock  propagation.  Chitty’s  work  infers  most 
of  the  Limestone  strength  and  compressibility  information  needed  to  construct 
models  such  as  SFT1  based  solely  on  porosity  and  RMR.  These  relations  are 
generic  for  limestone,  but  we  note  that  MQ  data  are  included  in  the  data  base 
used  to  generate  the  generic  relations.  We  wrote  an  automated  procedure  to 


4  Vaughan,  ‘Flex  User’s  Manual,  Version  1-J9.I,’,  Weidlinger  Associates,  2005. 

5  Mould,  J.  C.  and  H.S.  Levine,  “  A  Rate-Dependent  Three  Invariant  Softening  Model  for  Concrete,” 
Mechanics  of  Materials  and  Structures,  Voyiadjis,  G.  et.  al.  ,eds,  Elsevier,  Amsterdam,  1994. 

6  Chitty,  D.E.,  “A  Unified  Methodology  for  Estimation  of  Rock  mass  Properties  and  Associated 
Uncertainties  to  Support  Characterization  of  Inaccessible  Sites,”  DTRA  Test  Operations  Technology  and 
Test  support,  June  2005. 
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generate  SFT1  model  parameters  based  on  the  Chitty  correlations,  augmented  by 

some  further  assumptions  on  rate  effects,  extension  to  compression  ratios,  etc. 

Finite-element  modeling  and  ground-shock  parameters 
The  finite  element  analyses  for  ground-shock  simulations  were  performed  using  the  8- 
noded  hexahedron  element  of  the  FLEX  code.  The  box  shaped  study  area  was 
discretized  into  3’x  3’x  2’  finite  elements.  For  the  region  near  the  explosion  source, 
however,  a  finer  finite  element  mesh  was  used  to  resolve  the  explosive  charge.  The 
explosive  source  and  adjacent  region  were  modeled  using  the  FUSE  capability7  of  the 
FLEX  code.  The  source  was  a  51.875  inch  diameter  sphere  of  Nitromethane  -  the  same 
as  in  the  DG- 1  and  DG-2  experiments.  This  was  modeled  using  the  JWL  equation  of  state 
for  Nitromethane.  In  general,  this  approach  has  proved  satisfactory  for  ground  shock 
calculations.  For  this  exercise,  the  source  was  placed  in  the  center  of  the  characterized 
region  of  the  test  bed  rather  than  at  its  perimeter.  The  top  face  of  the  model  was  free, 
while  the  other  five  faces  were  restrained  using  absorbing  boundary  conditions.  Each 
analysis  was  performed  for  a  period  of  10  milliseconds  and  the  peak  hydrostatic  pressure, 
peak  compressive  principal  stress  and  peak  velocity  magnitude  during  this  period  were 
recorded  for  each  element.  These  aspects  of  the  calculation  remained  unchanged  for  all 
realizations  of  the  site 

Ground-shock  realizations  -  We  have  performed  two  series  of  stochastic  finite- 
element  simulations  using  the  WAI  FLEX  code  and  the  SFT1  constitutive  model,  and  the 
results  are  summarized  herein.  As  a  first  try,  we  generated  20  realizations  of  the  site 
based  on  the  vmodel  as  described  in  a  previous  section.  Each  realization  was  then  used  in 
a  ground-shock  simulation  calculation  subjected  to  the  same  energy  source.  Important 
free-field  measures,  such  as  peak  pressure,  peak  principal  stress,  peak  velocity,  etc.  were 
monitored  throughout  the  site.  Hence,  we  obtained  20  ground-shock  realizations  (free- 
field  environments)  that  corresponded  to  the  site  realizations.  We  denote  the  results  from 
this  series  as  “WAI  site”,  for  lack  of  a  better  term. 

In  the  second  series,  we  repeated  the  exercise  for  a  slightly  different  stochastic  site  model 
(derived  based  on  a  slightly  different  interpretation  of  the  same  database;  see  Sandia 
report).  Again,  the  same  finite-element  set  up  and  energy  source  was  used  so  that  the  only 
variable  was  the  site  property  realization.  These  results  are  referred  to  as  the  “S2  site”. 

We  have  also  performed  a  reference  calculation:  The  free-field  environment  was 
predicted  based  on  a  deterministic  site  model  which  corresponds  to  the  uniform,  layered 
site  geometry  used  by  the  community  in  the  DG2  exercise.  We  call  this  reference  the 
“layered  site” 

To  give  a  flavor  of  what  a  stochastic  ground-shock  environment  is  like,  we  compare  in 
Figure  16  one  realization  for  the  WAI  site  (sub-figure  (a))  with  one  for  the  S2  site  (sub- 


7  Sandler,  I.S.,  “A  new  Computational  procedure  for  wave  propagation  problems  and  a  new  procedure  for 
non-reflecting  boundaries,”  Computer  Methods  in  applied  mechanics  and  engineering,  Vol.  164,  1998. 
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figure  (b)).  Sub-figure  (c)  shows  the  environment  from  the  (deterministic)  layered  site, 
introduced  for  reference.  The  attribute  plotted  is  the  peak  pressure,  in  MPa.  For  clarity, 
the  figure  features  a  center-cut  of  the  volume  (viz.,  an  Easting  slice)  and  two  background 
cuts  for  3D  effects. 

The  peak  pressure  is  highest  at  the  source,  as  expected,  which  is  indicated  by  the  red 

O 

square  at  the  center  of  each  sub-figure.  But  since  the  source  pressure  is  high  (~  3,500 
Mpa),  it  dominates  the  rest  of  the  field  so  that  in  this,  the  default  magnitude  scale, 
differences  among  the  sub-figures  cannot  be  discerned  readily. 


Figure  16.  Comparing  a  ground-shock  realization  for  the  WA1  site  and  the  S2  site;  the 
result  from  the  deterministic  layered  site  serves  as  reference  purpose  only. 


8  Note  that  while  the  source  is  modeled  as  a  cube  (by  the  FUSE  module),  the  ground-shock  environment 
further  away  from  the  source  is  roughly  spherically  symmetrical  as  one  would  expect,  despite  the  presence 
of  layers. 
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To  see  the  lower-pressure  regions  more  clearly,  we  limit  the  magnitude  to  100  MPa,  and 
the  comparison  is  re-plotted  as  Figure  17.  The  layered  site  in  sub-figure  (c)  exhibits  a 
fairly  spherically  symmetric  environment,  which  is  only  tempered  by  the  presence  of 
layers.  On  the  other  hand,  the  environments  at  the  stochastic  sites  are  much  less 
symmetrical.  There  is  a  strong  angular  dependence  in  the  featured  cut  -  shades  of  a 
shooting  star,  versus  a  fireball.  The  “weakly  connected”  layered  features  in  the 
background  cuts  become  more  prominent,  reminiscent  of  the  layered  features  in  the 
stochastic  site-property  fields. 

To  see  yet  more  details  of  the  lower-pressure  regions,  the  cut-off  is  lowered  to  25  MPa 
and  the  results  presented  in  Figure  18.  The  fireball  becomes  more  ragged,  and  the  layers 
“connect”,  become  stronger  in  presence,  and  show  more  discontinuity  especially  when 
compared  with  the  reference.  The  WAI  site  and  S2  site  results  are  fairly  similar  in  general 
features,  albeit  they  can  be  very  different  when  compared  location  by  location  -  hence, 
stochasticity.  At  a  particular  location,  the  variation  from  realization  to  realization  can  be 
large  and  difficult  to  grasp,  exactly  the  experience  we  had  with  site  properties. 


23 


* 

Lower  value 

I  ^ 

en 


* 


(a)  WAI  stochastic  site,  realization  1 


(b)  S2  stochastic  site,  realization  1 


Figure  17.  Same  as  Figure  16,  but  with  the  magnitude  scale  limited  to  100  MPa. 


**•»»**»  r«*wfW 

(b)  S2  stochastic  site,  realization  1 


Figure  18.  Same  as  Figures  16  or  17,  but  with  the  magnitude  scaled  limited  to  25  MPa. 
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Figure  19.  Sample  mean  based  on  realizations  for  the  WAI  site.  Attribute  plotted  is  peak 

pressure,  in  MPa. 

Quantify  Local  and  Spatial  Variability  of  Ground-Shock  Response 


The  purpose  in  generating  the  sets  of  alternative  realizations  is  so  we  can  extract 
information  from  them  to  quantify  the  uncertainty  in  the  ground-shock  environment. 
Qualitatively,  one  can  get  a  feel  of  the  uncertainty  by  showing  the  realizations  visually 
one  at  a  time  but  in  rapid  succession,  like  frames  of  an  animated  cartoon9.  Quantitatively, 


9  Like  an  animated  cartoon,  successive  realizations  must  be  similar  enough  to  allow  the  eye  to  catch 
gradual  changes.  The  animated  display  of  realizations  allows  one  to  distinguish  areas  that  remain  stable 
over  all  realizations  (low  uncertainty)  from  those  with  large  fluctuations  between  realizations  (high 
uncertainty). 
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(a)  Same  as  previous  chart,  but  focusing  on  lower  magnitudes 


Figure  20.  Same  as  Figure  19  but  with  the  magnitude  scale  limited  to  100  MPa  to  show 

details  of  the  lower-pressure  regions. 

there  are  two  aspects  to  the  uncertainty:  Local  (point)  and  spatial  (global).  We  discuss 
them  further  in  the  following. 
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Local  uncertainty  is  quantified  by  the  sample  mean  and  variance,  which  are  typically 
combined  to  derive  a  (Gaussian-type)  confidence  interval  centered  on  the  mean  value10. 
With  more  assumptions,  the  sample  mean  and  variance  can  be  parlayed  into  a  local 
probability  distribution. 

Figure  19  gives  the  sample  mean  of  the  ground-shock  realizations  for  the  WAI  site. 
Again,  because  of  the  dominance  of  the  source  region  on  peak  value,  one  only  sees  the 
source  in  this  magnitude  scale  and  not  much  else.  We  circumvent  this  difficulty  as  before 
by  limiting  the  scale  to  100  MPa,  and  the  results  are  shown  in  Figure  20.  In  the  lower 
pressure  region,  the  sample  average  is  fireball  like  -  much  like  its  deterministic  (layered- 
site)  counterpart  given  previously  in  Figure  17(c). 


Figure  21.  Sample  standard  deviation,  WAI  site,  default  magnitude  scale.  Source  region 
has  zero  deviation  since  it  is  deterministic  -  same  source  is  used  for  all  realizations. 


10  For  example,  the  95%  confidence  interval  is  taken  as 

Pr-{Z(w)  E  [z  *  ( u )  -  2 o(u),  z  *  (u)  +  2 cr(w)]}=  0.95 

where  z*(u)  is  the  sample  mean  and  O'2  (u)  is  the  sample  variance  at  location  u  computed  from  the 
stochastic  realizations 
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Likewise,  we  show  the  sample  standard  deviation  in  Figures  21  and  22  in  two  magnitude 
scales.  The  plot  in  default  scale,  Figure  21,  is  dominated  by  the  near-source  region  - 
which  has  the  largest  standard  deviation.  Note  that  the  source  itself  has  zero  standard 
deviation  since  it  is  deterministic!  In  the  lower  magnitude  scale,  Figure  22,  the 
heterogeneity  shows  up  very  strongly. 


Figure  22.  Same  as  Figure  21,  but  limited  to  values  25  MPa  and  less. 

The  WA1  site  results  are  further  compared  with  the  S2  site  results  in  Figure  23,  to  give  an 
idea  of  the  sensitivity  to  slight  variations  in  the  site  vmodel.  A  couple  of  points  are  noted: 

(1)  Higher  pressure  regions  (other  than  the  source  region  itself,  which  is  deterministic) 
have  higher  standard  deviation 

(2)  Spatial  distributions  of  the  two  stochastic  sets  are  similar;  differences  are  reminiscent 
of  those  observed  in  a  comparison  of  the  site  realizations  -  S2  has  a  stronger  layered 
structure. 
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Hence,  it  is  desirable  to  use  a  measure  for  variability  that  is  independent  of  the  ground- 
shock  magnitude.  This  prompts  us  to  look  to  the  coefficient  of  variation,  which  is  the 
sample  standard  deviation  normalized  by  the  sample  mean  and  commonly  denoted  as 
COV.  As  can  be  seen  from  the  comparison  in  Figure  24,  the  COV  runs  between  0  and  1, 
with  higher  COVs  distributed  near  the  source  and  lower  COVs  further  away  from  the 
source.  The  S2-site  COVs  exhibit  a  stronger  layer  structure;  the  WAI-site  results  are 
more  salt-and-pepper.  The  histograms  on  the  right  side  of  the  figure  indicate  that  the  two 
COV  distributions  are  in  fact  very  similar.  Mean  value  of  the  WAI-site  COV  is  slightly 
larger,  indicating  slightly  greater  variability  in  the  WAI-site  ground-shock  environment. 
This  is  a  direct  result  of  the  higher  variability  in  the  site  model,  as  was  mentioned  in  a 
previous  section. 

Finally,  using  the  sample  mean  and  standard  deviation  (or  COV),  bounding  fields  such  as 
those  based  on  mean  ±  n*standard  deviation,  where  n  is  an  integer,  can  be  generated.  For 
example,  Figure  25  gives  the  (mean  ±  l  standard  deviation)  bounds.  The  uncertainty  in 
ground-shock  at  any  location  in  the  field,  in  this  case  the  peak  pressure,  is  thus 
delineated. 

It  is  obvious  that  other  ground-shock  attributes  can  be  processed  in  identical  fashion.  This 
is  not  done  due  to  time  and  budget  limitations. 

Spatial  uncertainty  of  ground-shock,  like  that  of  the  site  properties,  requires  a  measure  of 
the  joint  uncertainty  about  ground-shock  values  at  several  locations  taken  together  -  e.g., 
the  probability  of  occurrence  of  a  string  of  high  peak  stresses  or  low  stresses.  Such 
information  can  be  extracted  from  the  stochastic  fields  using  the  same  geostatistics 
techniques  used  for  site  characterization.  Alternately,  the  uncertainty  information  can  be 
summarized  using  different  types  of  metrics  including:  Probability  maps,  quantile  maps, 
and  maps  of  spread.  These  directions  will  be  explored  in  a  subsequent  study. 
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(b)  S2  stochastic  site 
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Figure  23.  Comparison  of  sample  standard  deviations,  realizations  for  WAI  site  and  the 

S2  site. 
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Figure  24.  Comparison  of  COVs  based  on  realizations  for  the  WAI  site  and  the  S2  site. 
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(a)  WAI  stochastic  site 


Figure  25.  Upper  and  lower  bounds  of  peak  pressure,  as  measured  by  (sample  mean  ±  1 
standard  deviation). 
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As  a  final  point,  figures  26  and  27  compare  velocity  waveforms  from  several  stochastic 
realizations  to  that  from  the  conventional  layered  approximation.  Note  the  smoother 
appearance  and  slightly  different  structure  of  the  heterogeneous  sites  relative  to  the 
layered  site.  Though  not  analyzed  in  this  report,  we  find  these  differences  very 
interesting. 


Figure  26.  Radial  velocities  13.5  ft.  below  charge. 
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Figure  27.  Radial  velocities  34.6  ft.  below  charge. 
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HYPERLINKED  TEXTS 


Appendix  A.  Selected  Sandia  Results 

We  compare  the  vmodel  developed  herein  with  those  discussed  in  the  Sandia  report,  for 
reference  purposes.  All  variograms  have  been  normalized  to  have  a  sill  value  of  L  for 
easy  comparison.  The  model  labeled  WAI  is  the  model  developed  from  data  from  the  33 
RPSS  soundings  as  described  in  the  main  text.  The  model  labeled  Sandial  (SI) 
corresponds  to  that  given  in  the  Sandia  report  based  on  23  RPSS  soundings  (Phase  1  site 
exploration.  Table  3  of  Sandia  report),  whereas  the  Sandia2  (S2)  model  corresponds  to 
that  given  in  Table  6  of  said  report  for  a  subset  of  the  volume  of  interest. 

With  reference  to  Figure  A-l,  the  vertical  variograms  for  WAI  and  SI  are  identical,  for 
all  practical  purposes.  The  horizontal  variograms  differ  in  two  important  aspects:  the 
range  and  asymptote.  In  particular,  the  effective  horizontal  range  is  about  100+  ft  for 
WAI  and  only  15  for  SI.  Hence,  one  expects  the  SI  field  to  be  choppier  horizontally. 
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(b)  Vertical  variogram 


Figure  A-l.  Comparison  of  vmodels  generated  herein  and  in  said  reference. 
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Compared  with  the  WAI  model,  the  S2  model  has  a  smaller  vertical  range  (8  ft  versus  20 
ft),  and  a  larger  horizontal  range  (over  200  ft,  versus  100+  ft).  Hence,  one  expects  the  S2 
field  to  be  more  choppy  vertically,  but  very  uniform  horizontally.  These  effects  are 
confirmed  by  the  corresponding  simulated  fields  given  in  Figure  A-2.  Note  the 
correspondence  between  large  effective  horizontal  ranges  (Sandia2)  and  more  uniform 
layering  structure  in  the  simulated  field. 
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Figure  A-2.  Comparison  of  realizations  based  on  the  WAI  vmodel  and  the  Sandial  and 

Sandia2  vmodels  given  in  Figure  A-l. 
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